## load necessary libraries
##
library(Cairo)
library("RSVGTipsDevice")

## initialize device
##
CairoWin(width = 7, height = 7, pointsize = 12, record = getOption("graphics.recorecord"),
rescale = c("R", "fit", "fixed"), xpinch, ypinch, bg = "transparent", canvas = "white", 
gamma = getOption("gamma"), xpos = NA, ypos = NA, buffered = getOption("windowsBuffered"),
restoreConsole = TRUE)
par(mfrow=c(2,1))
## pie(rep(1,32), col=heat.colors(32))
colors<-heat.colors(32)

## BACK TO BUSINESS
## first of all - load the data
##
sample <- read.csv(file="senin.devtime.CSV", header=TRUE, na.strings = "null", sep=",", quote="\"'", dec=".")
# make vectors of dates and values for easier analysis / plotting
sample.t  <- as.Date(sample[,1])
sample.v <- as.vector(sample[1-dim(sample)[1],9])
# cure n/A's with 0
for(i in 1:length(sample.v)){ if(is.na(sample.v[i])) { sample.v[i] = 0} }

## plot the data in simple manner, just a plot
##
plot(tsdata.values, type="l", lwd=2, col=colors[14], xaxt = "n", ylab="devTime, hrs", main="Development Time")
abline(h=mean(tsdata.values), lwd=6, col=paste(substr(colors[15],1,7),"88",sep=""))
#barplot(data.values, rep(1, length(data.values)), width=5, col=colors[2], xaxt = "n", ylab="devTime, hrs", main="Development Time")
# extract month ticks and make labels
tmp <- format(tsdata.times, "%b %y")
month.labels <- unique(format(tsdata.times, "%b %y"))
month.at <- rep(-1,length(month.labels))
j <- 1
for(i in 1:length(tmp)){
 if( (month.labels[j] == tmp[i]) && (month.at[j] == -1 ) ){ month.at[j] = i; j <- j+1; }
}
axis(1, at=month.at, labels=month.labels, tcl=-1)
# make days axis
days.at <- c(1:length(tsdata.times))
days.labels <- format(tsdata.times, "%d")
#
axis(1, at=days.at, labels=FALSE, tcl=-0.1)



## FFT stuff goes in here
##
tsdata.fft <- fft(tsdata.values)
tsdata.f.Nyquist <- 1 / 2 / sampling.delta
tsdata.f <- tsdata.f.Nyquist * c(seq(length(tsdata.times)/2), -rev(seq(length(tsdata.times)/2))) / (length(tsdata.times)/2)
## plot it
plot(tsdata.f[2:length(tsdata.f)], Mod(tsdata.fft)[1:length(Mod(tsdata.fft))]/length(Mod(tsdata.fft)), xlim=c(-2,2), type='o', lwd=2, 
 col=heat.colors(32)[15], xlab="Frequency, Hertz", ylab="Power", main="Simple spectral analysis")














## plot the development time data in simple manner, just a plot
##
hist(diff(tsdata.values),prob=T,ylim=c(0,0.5),xlim=c(-6,6),col=colors[14],
 main="Development Time distribution properies")
lines(density(diff(tsdata.values)),lwd=3, col=colors[7])

## plot the development time QQ normal plot
##
qqnorm(diff(tsdata.values), col=colors[14], main="Normal Q-Q Plot of Development Time")
abline(0,1,col=colors[7], lwd=2)


## Find the sample period:
#delta <- ecg$t[2] - ecg$t[1]
#
# TODO: I think this must be fixed, since we are getting a day as the period
# 
delta <- as.numeric(tsdata.times[2]-tsdata.times[1])


## Calculate the Nyquist frequency:
f.Nyquist <- 1 / 2 / delta

## Calculate the frequencies.  (Since ecg$t is in seconds, delta
## is in seconds, f.Nyquist is in Hz and ecg$freq is in Hz)
## (Note: I may be off by 1 in indexing here ????)
data.freq <- f.Nyquist*c(seq(length(data.values)/2), -rev(seq(length(data.values)/2)))/(length(data.values)/2)
plot(data.freq, type='l', main='FFT of ECG vs frequency', xlab='Frequency [Hz]') 

## Plot fft vs frequency
plot(Mod(fft) ~ freq, data=ecg, type='l', log='y', main='FFT of ECG vs frequency', xlab='Frequency [Hz]') 




#################################################
## component analysis
##
## library(ts)
## plot(data.values,type="l")
## tui.1 <- filter(data.values,filter=rep(1/5,5))
## tui.2 <- filter(data.values,filter=rep(1/25,25))
## tui.3 <- filter(data.values,filter=rep(1/81,81))
## lines(tui.1,col="red")
## lines(tui.2,col="purple")
## lines(tui.3,col="blue")


#par(mfrow=c(2,2))
#spectrum(data.ts)
#spectrum(data.ts, spans=3)
#spectrum(data.ts, spans=c(3,3))
#spectrum(data.ts, spans=c(3,5))



## now plot the data
##
# extract ticks text
# make TS (time series) object from the data
data.ts <- ts(ts(data.values),start=c(data[2,4], data[2,8]), freq=7)
plot(data.ts, xaxt="n")
axis.Date(1, at=seq(as.Date(data[2,1]), as.Date(data[x_ticks.num,1]), "weeks"))
axis.Date(1, at=seq(as.Date(data[2,1]), as.Date(data[x_ticks.num,1]), "days"), labels = FALSE, tcl = -0.2)
print(data.ts)
